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Abstract 

We consider two-component fermions on the lattice in the unitarity limit. This is an idealized 
limit of attractive fermions where the range of the interaction is zero and the scattering length 
is infinite. Using Euclidean time projection, we compute the ground state energy using four 
computationally different but physically identical auxiliary-field methods. The best performance 
is obtained using a bounded continuous auxiliary field and a non-local updating algorithm called 
hybrid Monte Carlo. With this method we calculate results for 10 and 14 fermions at lattice 
volumes 4 3 , 5 3 , 6 3 , 7 3 , 8 3 and extrapolate to the continuum limit. For 10 fermions in a periodic cube, 
the ground state energy is 0.292(12) times the ground state energy for non-interacting fermions. 
For 14 fermions the ratio is 0.329(5). 
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I. INTRODUCTION 



The unitarity limit describes attractive two-component fermions in an idealized limit 
where the range of the interaction is zero and the scattering length is infinite. The name 
refers to the fact that the S-wave cross section saturates the limit imposed by unitarity, 
0o(&) — 4:71 /k 2 , for low momenta k. In the unitarity limit details about the microscopic 
interaction are lost, and the system displays universal properties. Throughout our discussion 
we refer to the two degenerate components as up and down spins, though the correspondence 
with actual spin is not necessary. At sufficiently low temperatures the spin-unpolarized 
system is believed to be superfluid with properties in between a Bardeen-Cooper-Schrieffer 
(BCS) fermionic superfluid at weak coupling and a Bose-Einstein condensate of dimers at 

ndri 

strong coupling [l|, [2j, [3]. 

In nuclear physics the phenomenology of the unitarity limit is relevant to cold dilute 
neutron matter. The scattering length for elastic neutron-neutron collisions is about —18.5 
fm while the range of the interaction is comparable to the Compton wavelength of the pion, 
m^ 1 pa 1.4 fm. Therefore the unitarity limit is approximately realized when the interparticle 
spacing is about 5 fm. While these conditions cannot be produced experimentally, neutrons 



at around this density may be present in the inner crust of neutron stars 

The unitarity limit has been experimentally realized with cold atomic Fermi gases. For 
alkali atoms the relevant length scale for the interatomic potential at long distances is the 
van der Waals length £ v dw- If the spacing between atoms is much larger than £ v dw, then 
to a good approximation the interatomic potential is equivalent to a zero range interaction. 
The scattering length can be manipulated using a magnetic Feshbach resonance 

Baa 
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121 ] . This technique involves a molecular bound state in a "closed" hyperfine 
channel crossing the scattering threshold of an "open" channel. The magnetic moments 
for the two channels are typically different, and so the crossing can be produced using an 
applied magnetic field. 

At zero temperature there are no dimensionful parameters other than particle density. 
For iVj up spins and N± down spins in a given volume we write the energy of the unitarity- 
limit ground state as N . For the same volume we call the energy of the free non- 
interacting ground state E^ he ^ , and the dimensionless ratio of the two £jv T) jvj> 



6v T ,JV; - E° N ^ j7V| / E^*^ . ( 1 ) 



The parameter £ is denned as the thermodynamic limit for the spin-unpolarized system, 



lim £ NN . 

N^oo 



(2) 



Several experiments have measured £ from the expansion of Li and K re 



larmonic trap. Some recent measured values for £ are 0.51(4) 12 1. 0.46 + 5 13], and 0.32^3 



eased from a 



ll| . The disagreement between these measurements and larger values for £ reported in 



earlier experiments |6J, O, Q suggest further work may be needed. 

There are numerous analytic calculations of £ using a variety of techniques such as BCS 
saddle point and variational approximations, Pade approximations, mean field theory, den- 
sity functional theory, exact renormalization gro up, dimensional e-expansions, and large- N 
expansions Q Q, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26]. The values for £ range from 



0.2 to 0.6. Fixed-node Green's function Monte Carlo simulations for a periodic cube have 



found £ NjN to be 0.44(1) for 5 < iV < 21 [27( and 0.42(1) for larger N [28 



291 ] . A restricted 



path integral Monte Carlo calculation finds similar results 30[], and a sign-restricted mean 

field lattice calculation yields 0.449(9) 3l| . 

There have also been simulations of two-component fermions on the lattice in the unitarity 
imit at nonzero temperature. When data are extrapolated to zero temperature the results of 
32 , 33| produce a value for £ similar to the fixed- node results. The same is true for 34 



35], 



; hough with significant error bars. The extrapolated zero temperature lattice results from 
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371 ] established a bound, 0.07 < £ < 0.42, while more recent lattice calculations yield 



0.261(12) 



38 



391 ] . The work of [38j, [39] includes a comparison of two lattice algorithms 



similar to the ones discussed in this paper. 

In 40(] £,n,n was calculated on the lattice using Euclidean time projection for iV = 
3, 5, 7, 9, 11 and lattice volumes 4 3 , 5 3 , 6 3 . From these small volumes it was estimated that 
£ = 0.25(3). More recent results using a technique called the symmetric heavy-light ansatz 
found similar values for £jv,iv at the same lattice volumes and estimated £ = 0.31(1) in the 
continuum and thermodynamic limits {41 ]. 

For lattice calculations of £n,n at fixed N, systematic errors due to nonzero lattice spacing 
can be removed by extrapolating to the continuum limit. At unitarity where the scattering 
length is set to infinity, this continuum extrapolation corresponds with increasing the lattice 
volume as measured in lattice units. In this paper we present new lattice results using 
the Euclidean time projection method introduced in [4(3] and extrapolate to the continuum 
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limit. 

In earlier work 0, singular matrices were encountered that prevented calculations at 
larger volumes. We discuss this and related issues by considering four different auxiliary-field 
methods which reproduce exactly the same Euclidean lattice model. These four methods 
include a Gaussian-integral auxiliary field, an exponentially-coupled auxiliary field, a discrete 
auxiliary field, and a bounded continuous auxiliary field. By far the best performance is 
obtained using the bounded continuous auxiliary field and hybrid Monte Carlo [43]. With 
this method we calculate the ground state energy for 10 and 14 fermions at lattice volumes 
4 3 , 5 3 , 6 3 , 7 3 , 8 3 and extrapolate to the continuum limit. 



II. EUCLIDEAN TIME LATTICE FORMALISM 



Throughout our discussion of the lattice formalism we use dimensionless parameters and 
operators corresponding with physical values multiplied by the appropriate power of the 
spatial lattice spacing a. In our notation the three-component integer vector n labels the 
lattice sites of a three-dimensional periodic lattice with dimensions L 3 . The spatial lattice 
unit vectors are denoted I — i, 2, 3. We use n t to label lattice steps in the temporal 
direction. L t denotes the total number of lattice time steps. The temporal lattice spacing 
is given by at, and at = at /a is the ratio of the temporal to spatial lattice spacing. We also 
define h = at/ (2m), where m is the fermion mass in lattice units. 

We discuss four different formulations of the Euclidean time lattice for interacting 
fermions: Grassmann path integrals with and without auxiliary field and transfer matrix 
operators with and without auxiliary field. For any spatial and temporal lattice spacings 
these four formulations produce exactly the same physics, as indicated in Fig. [IJ We follow 
the order of discussion indicated by the numbered arrows in Fig. [U 



A. Grassmann path integral without auxiliary field 

Let Q and c* be anticommuting Grassmann fields for spin i. The Grassmann fields are 
periodic with respect to the spatial lengths of the L 3 lattice, 

q(n + LI, n t ) = Ci(n + L2, n t ) = q(n + L3, n t ) = Ci(n, n t ), (3) 
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Grassmann path integral 
without auxiliary field 

exp [-S{c c*)] 



#1 



#2 



#3 



Grassmann path integral 
with auxiliary field 

exp [— S(c, c*, s)) 



(four versions) 



Transfer matrix operator 
without auxiliary field 

M 



Transfer matrix operator 
with auxiliary field 

M{s, n t ) 



(four versions) 



#4 



Monte Carlo calculations 

FIG. 1: (Color online) Schematic diagram showing four equivalent Euclidean time lattice formula- 
tions. The numbered arrows indicate the order of discussion in the text. 



and antiperiodic along the temporal direction, 

d(n, n t + L t ) = -Ci(n, n t ). 
We write DcDc* as shorthand for the integral measure, 

DcDc* = J l dci(n, n t )dc*(n, n t ) 

n,nt,i=1,l 

The Grassmann spin densities, pj and p±, are defined as 

p T (n,n 4 ) = c*(n,n t )ci(n,nt), 

Pl(n,n t ) = cj(n,n 4 )q(n,n 4 ). 
We consider the Grassmann path integral 



Z = J DcDc* exp [-S (c, c*)] , 
S(c,c*) = S iree (c,c*) + Ca t 22pi(n, n t )pi(n, n t ) 



(4) 



(5) 



(6) 
(7) 

(8) 
(9) 



The action S(c,c*) consists of the free fermion action 



Sf ree (c, c*) = ^ [c* (n, n t )cj(n, n t + 1) - (1 - 6h)c*(n, n t )d(n, n t )} 

n,rat,i=t,J, 

-fr ^ k*(™> n *) c i(™ + ^> n *) + c*(n,n t )ci(n- l,n t ) 

n,nt,«=T,J. 2=1,2,3 



(10) 



and an attractive contact interaction between up and down spins. We take this action as 
the definition of the lattice model. All other formulations introduced later will be shown 
to be exactly equivalent. 

We use a fermion mass of 939 MeV and lattice spacings a = (50 MeV) -1 , a t = (24 
MeV) -1 . In dimensionless lattice units the corresponding parameters are m = 18.78 and 



at = 2.0833. These are the same values used in 
scales for dilute neutron matter. 



40] and is motived by the relevant physical 



B. Interaction coefficient C for general spatial and temporal lattice spacings 



We use Liischer's formula 



44 



45 



1461 . 147j to relate the coefficient C to the S-wave scatter- 
ing length. We consider one up-spin particle and one down-spin particle in a periodic cube 
of length L. Liischer's formula relates the two-particle energy levels in the center-of-mass 



frame to the S-wave phase shift, 

J9COt5 (p) 



7lL 



Lp 
2^ 



where S(r]) is the three-dimensional zeta function, 



S(V) 



lim 

A^oo 



^ 6(A 2 - n 2 ) 

22 — ~ 47rA 



For small momenta the effective range expansion gives 

pcot^oGo) P 



1 1 2 
h 7; r oP + 

Oscatt * 



where a scat t is the scattering length and r is the effective range. 
In terms of r], the energy of the two-particle scattering state is 



E- 



pole 



p 



7] (Itx 
m\L 



(11) 



(12) 



(13) 



(14) 
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oo 

= z >o-o<; 



FIG. 2: Sum of bubble diagrams contributing to two-particle scattering. 

We compute the location of the two-particle scattering pole by summing the bubble diagrams 
shown in Fig. [2J After writing the sum of bubble diagrams as a geometric series, the relation 
between C and E^ Q \ e becomes 

- -— = lim — > - - , 15 

Ca t l^ooL 3 ^ e -E polc a t -i + 2a t uj(2Tik/L) - a 2 u 2 (27ik/L) 

where 

U (P) = — Y] (1-cospi). (16) 
m L — 4 

/=1,2,3 

These expressions were first derived in [481 ] . In the unitarity limit the scattering length a sca tt 
is set to infinity and all other effective range parameters are set to zero. From Eq. f[T3|) 
and ( TTTT) this is equivalent to setting S (77) = 0. To minimize lattice discretization errors we 
pick the root of S (if) closest to zero, rj = —0.095901. After selecting some large value for 
L, we use Eq. (1141) and (fT5l) to determine the unitarity limit value for the coefficient C. For 
our chosen lattice parameters C is —7.4415 x 10 -5 MeV -2 in physical units or —0.18604 in 
dimensionless lattice units. 



C. Transfer matrix operator without auxiliary field 



We convert the Grassmann path integral into the trace of a product of transfer matrix 
operators using the exact correspondence Q, Q 



Tr : F if , 



aUn'),aAn] 



al(n'),ai(n) : j 



/ DcDc* exp <j ^ ^2 c **^' n ^ t Ci ^' n ^ ~ Ci ^' nt + 1 )1 f 

I n t =0 n,i J 



Lt-1 



X 



II [c*,(n',n t ),Ci(n } n t )} 



(17) 



n t =0 
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where Ci(n, L t ) = — Cj(n, 0). We use Oj(n) and a|(n) to denote fermion annihilation and 
creation operators respectively for spin i at lattice site n. The : symbols in (TTTjl indicate 
normal ordering. Let us define 
3 



//iv. .. = — ^ at(n)ai(n) - ^— ^ ^ [al(n)ai(n + I) + a\(n)ai(n - I) (18) 

n,j=T,J, fM=t,J. 2=1,2,3 



and 



pf a (n) = o}(n)a T (ft), 



Then 

Z = Tr (M Lt ) , 
where M is the normal-ordered transfer matrix operator 



M =: exp 



n 



(19) 
(20) 

(21) 
(22) 



D. Grassmann path integral with auxiliary field 



We consider four different auxiliary-field transformations labelled with the subscript j = 
1, 2, 3, 4. Let us define the Grassmann actions 

Sj (c,c*,s) = S ivcc (c,c*) -^Aj [s(n,n t )) ■ [p^(n,n t ) + ^(n,n t )] , (23) 
and Grassmann path integrals 



n,n t 



£j = Y\_ J djs(n,n t ) J DcDc* exp [— Sj (c, c*, s)} 

n,n t 



(24) 



where s is a real- valued function over all lattice sites. For j = 1, we use a Gaussian-integral 
transformation similar to the original Hubbard-Stratonovich transformation 5l|, [52] , 

i p+oo 



dis(n, n t ) 



^7T 



ds(n,n t )e 2 



s 2 (n,n t ) 



A\ [s(n,n t )] = ^-C 1 a t s{n,n t ). 



For j = 2 we use an auxiliary field with exponential coupling to the fermion densities, 

i r+oo 

—= / ds{n,n t )e-^ s2 ^ nt \ 

V27T J-oo 



d 2 s(n, n t ) 



(25) 
(26) 

(27) 
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A 2 [s{n,n t )\ = (1 - 6/1) 



iV-C20its(ri,nt)+ ° 2 " t 



This is the auxiliary- field transformation used in 36|, [37, 



40 



42 



(28) 



481 ]. For j = 3 we use a 



discrete auxiliary field taking values —1 or +1 at each lattice site, 

J d 3 s(n,n t ) = ~ ^ 



s(ra,n t )=±l 

A 3 [s(n,n t )] = yJ-C z a t s{n, n t ). 



(29) 



(30) 



The discrete Hubbard-Stratonovich transformation was first introduced by Hirsch in [53J 



and has been used in a number of recent lattice simulations 



32 



33, 



38 



39. 



5J] . For j = 4 



we introduce a new transformation which uses a bounded but continuous auxiliary field, 



d A s(n, n t ) 



1 

2^ 



ds(n, n t ) 



A 4 [s(n, nt)] = \/-C 4 a t sin [s(n, n t )\ . 



(31) 



(32) 



The motivation for the new transformation is discussed later when comparing and analyzing 
computational performance. 

For each j = 1, 2, 3, 4 the integral transformations are defined so that 



djs(n, n t )l = 1, 



djs(n, n t ) Aj [s(n, n t )] = 0. 



(33) 



(34) 



Since all even products of Grassmann variables commute, we can factor out the term in Zj 
involving the auxiliary field s at n, n t . To shorten the notation we temporarily omit writing 
n, nt explicitly. For each j = 1, 2, 3, 4 we find 

J djsexp [Aj (s) (p T + p x )] = J djs [l + Aj (s) (p T + pj + A) (s) p ]Pl ] 

= l + y djsA) (s)ptP| = exp J d j sA 2 j{s)p^p l . (35) 



For the four cases the integral of A 2 - is 



J disA\ (s) = -C x a u 
j d 2 s A 2 2 (s) = (1 - 6h) 2 (e- c '^-l), 



(36) 
(37) 



TABLE I: Filling sequence of momentum states for each spin. 



N 


additional momenta filled 


1 


(0,0,0) 


3 


(f>°>°>> <-x>°>°> 


5 


(0,^,0), (0,-^,0) 


7 


(0,0,f),<0,0,-f) 



J d 3 sA 2 3 0) = -C 3 a t , 



d±s A\{s) = —-C^a t . 



Therefore 

Z = Z\ = Z2 = -Zg = Zt± 

when the interaction coefficients satisfy 

(l-6h) 2 (e- c * a * - 1) 
C = d = ^ '- = C 3 

a t 



\c. 



(38) 
(39) 

(40) 
(41) 



III. TRANSFER MATRIX OPERATOR WITH AUXILIARY FIELD 



Using Eq. f[T7l) and fl40|) we can write Z as a product of transfer matrix operators which 
depend on the auxiliary fields, 



n,nt 



djs(n, n t ) 



TriMj^Lt-l) 



Mj{s,0)} 



(42) 



where 



Mj(s, n t ) =: exp \ -H ircc a t + ^ Aj [s(n, n t )} ■ \pf a (n) + pf a (n) 



(43) 



For the numerical lattice calculations presented in this paper we use the auxiliary-field 
transfer matrix operators in Eq. (|43j) . 



Let 



0,free\ 
N,N I 



be the normalized Slater-determinant ground state on the lattice for a non- 



interacting system of N up spins and N down spins. For each spin we fill momentum states 
in the order shown in Table [B We construct the Euclidean time projection amplitude 



Z N ,N{t) = \\ J djs(n, n t ) 

n,n t 



0,free 
N,N 



MjfaLt-l) 



0,free 
N,N I ' 



(44) 
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where t = L t a t . The result upon integration gives the same Z^^it) for each j = 1, 2, 3, 4. 

As a result of normal ordering, Mj(s,nt) consists of only single-particle operators inter- 
acting with the background auxiliary field and no direct interactions between particles. We 
find 

Mj(s,L t - 1) Mj-foO) Ufffi) = [detM,(M)] 2 , (45) 



,T,0,free 



where 

[Mi(s, t)] k , k = (pvlMjis, L t -l) Mj{s, 0) \p k ) , (46) 

for matrix indices k, k' — 1, • - • , N. \pk) , \vk') are single-particle momentum states compris- 
ing the Slater-determinant initial/final state. The single-particle interactions in Mj(s,n t ) 
are the same for both up and down spins, and this is why the determinant in Eq. (145!) is 
squared. Since the matrix is real-valued, the square of the determinant is nonnegative and 
there is no problem with oscillating signs. 



For the interacting lattice system at unitarity, we label the energy eigenstates I^/^tv) 

<k 

N.N 



with energies E% N in order of increasing energy. 



E N)N < Ejf iN < ■ ■ < E k N ^ N < ■ ■ ■ . (47) 

In the transfer matrix formalism these energies are defined in terms of the logarithm of the 
transfer matrix eigenvalue, 

M\^ N ) = e- E ^ at \^ k NtN ). (48) 
Let us define c k N N as the inner product with the initial free fermion ground state, 

4,^ = <^|^5 e )- (49) 

In the following we assume that c° NN , the overlap between free fermion ground state and 
the interacting ground state, is nonzero. 

Let us define a transient energy expectation value that depends on the Euclidean time t, 

tp u\ 1 i Z N N (t - a t ) , , 

E N , N {t) = — In — -— . (50) 

The spectral decomposition of Z^^{t) gives 

;w*) = El c ^v|V^-\ (si) 
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and at large Euclidean time t significant contributions come from only low energy eigenstates. 
For large t we find 



E N , N (t) - & Nl s + E 7 ' ' " e-W*-*M*. (52) 



For low energy excitations E^ N — E° N N is much smaller than the energy cutoff scale a t 1 
imposed by the temporal lattice spacing. Therefore 



\ C N,N 


2 

e 


, E N,N- E N,N) a t — I 


\r° 

\ N,N 


2 





E N , N (t) - K,n + EfrT ^ - K ^ (53) 

fc^O \ C N,N | 

The ground state energy N is given by the limit 

E° N>N = lim E NtN (t). (54) 



t— >oo 

IV. IMPORTANCE SAMPLING 

We calculate the ratio 

Zn,n\P ~ a t) / rrN 

using Markov chain Monte Carlo. For t = L t at, configurations for the auxiliary field are 
sampled according to the weight 

exp{-£/ i (a) + 21n[|detM i (s,L t a t )|]}, (56) 

where 

UUs) = Uo(s) = 

n,nt 

and 



Ul{s) = U2{s) = ^J2[s(n,n t )} 2 , (57) 



U 3 {s) = U 4 (s) = 0. (58) 



For j = 1, 2, 4 importance sampling is implemented using hybrid Monte Carlo 43]. This 
involves computing molecular dynamics trajectories for 

H j (s, P ) = ^\p(n,n t )] 2 + V j (s), (59) 

n,nt 

where p(n,nt) is the conjugate momentum for s(n,nt) and 

Vj(s) = E7 J -(s)-21n[|detM i (s,L t a t )l]- ( 60 ) 
The steps of the algorithm are as follows. 
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Step 1: Select an arbitrary initial configuration s°. 



Step 2: Select a configuration p° according to the Gaussian random distribution 



P [p°(n,n t )] oc exp <j ~ [p°(n,n t )\ 



Step 3: For each n,n t let 



p°(n,n t ) = p°(n,n t ) - 



^stcp 



ds(n, n t ) 



for some small positive e stcp . The derivative of Vj is computed using 

dV^s) _ dU 3 (s) 



ds(n,n t ) ds(n,n t ) det Mj ^ d [Mj] kl ds(n, n t ) 



2 ^ d&etMj <9[M 



(9s(n, n t ) 



Step 4: For steps i = 0, 1, iV ste p — 1, let 



dUAs) r.x n 9 [ M 

- 2 E[ M 7] ' ' 



-3\kl 



k,l 



' lk ds(n, n t )' 



s l+L (n,n t ) = s l (n,n t ) + e stcp p l (n, n t ), 
p l+1 (n, n t ) = p l (n, n t ) - e stcp 



dV 3 {s) 



ds(n, n t ) 



for each n, n t . 



Step 5: For each n,n t let 



p N «*>(n, nt) = p Nst ^(n, n t ) + ^ 

2 



Step 6: Select a random number r G [0, 1). If 



dV(s) 
ds(n, nt) 



r < exp [-H(s Nate »,p Nstep ) + H(s ,p )] 



(61) 



(62) 



(63) 



(64) 
(65) 



(66) 



(67) 



then set s° = s Nstcp . Otherwise leave s° as is. In either case go back to Step 2 to 
start a new trajectory. 

For j — 3, hybrid Monte Carlo cannot be used since the auxiliary field has only discrete 
values —1 or +1. In this case we use a local Metropolis accept/reject update. For each 
sweep through the entire lattice a small fraction of random lattice sites are selected. A new 
configuration s' is produced by flipping the sign of s at these sites. A random number r G 



13 



[0, 1) is selected and the new configuration s' is accepted or rejected based on the Metropolis 
condition, 

exp {-U 3 (s') + 2 In [|det M 3 (s', L t a t )\}} 
T ex P {-U 3 (s) + 2ln[\detM 3 (s,L t a t )\}} ' [ ' 

For all cases j = 1, 2, 3, 4 the observable we calculate for each configuration is 

[detM^s.L.a,)] 2 
By taking the ensemble average of Oj(s, L t at) we obtain 

ZN f~f (70) 

for t = L t a t . 



V. PRECISION TESTS AND PERFORMANCE COMPARISONS 
A. Precision tests for two particles 

We perform precision tests of the four auxiliary-field methods using the two-particle 
system at unitarity. For the system with one up spin and one down spin, we compute 
the dimensionless combination mL 2 Ei^(t) for L = 4 and L t = 6, 12. As explained in the 
previous section we use hybrid Monte Carlo for j = 1, 2, 4 and the Metropolis algorithm for 
j = 3. The simulations are performed using 16 processors starting with different random 
number seeds producing independent configurations. The final result is computed from the 
average of the individual processor results, and the standard deviation is used to estimate the 
error of the average. The linear space of states for the two-particle system is small enough 
that the transfer matrix without auxiliary field in Eq. (1221) can be used to calculate exact 
results for comparison. The results are shown in Table [III For each case, j = 1, 2, 3, 4, the 
auxiliary-field Monte Carlo results using Mj reproduce the exact results up to the estimated 
stochastic errors. 



B. Performance comparisons 

In our discussion E^ 1 ^ refers to the ground state energy for N up-spin and N down-spin 
free fermions on the lattice. More explicitly is the energy exponent of the free-particle 
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TABLE II: Comparison of exact transfer matrix results and j 
Carlo results for mL 2 Ei < i(t) with L = 4 and Lt = 6, 12. 



1,2,3,4 auxiliary- field Monte 



L t 


Exact result 


Mi 


M 2 


M 3 


M 4 


6 


-2.087 


-2.06(3) 


-2.10(3) 


-2.05(4) 


-2.10(3) 


12 


-2.817 


-2.77(3) 


-2.81(4) 


-2.80(4) 


-2.88(4) 



transfer matrix, 



: exp (-H hee a t ) : 



0,frcc 
JV,AT 



= exp 



/ i-iO.frcc \ 
[~ E N,N a t ) 



0,free 
N,N 



(71) 



for the same lattice parameters used to calculate £jv,jv(0- I n contrast with this lattice energy 



definition for E^f 1 ^, we define the Fermi energy .Ep purely in terms of particle density. For 



iV up spins and iV down spins in a periodic cube we let 



67T 



2/3 



~ 7.596 



2m 2m V, L 3 J 
For > 1 it is useful to define the dimensionless function 

EN,N{t) 



AT2/3 

mL 2 



£,N,N(t) 



E 



i0,free 
N,N 



(72) 



(73) 



Scale invariance at unitarity requires that in the continuum limit £,N,N{t) depends only on 
the dimensionless combination -^i- For fixed the Fermi energy Ep is proportional to 
and so we can regard ^N,N(t) as a universal function of E F t. This universal behavior 
provides a simple but nontrivial check of scale invariance in our lattice calculations. 

As one part of our analysis of computational performance we monitor the rejection rate 
at Step 6 of the hybrid Monte Carlo algorithm for j = 1,2, 4 as well as the Metropolis 
update for j = 3. The average likelihood of rejection is recorded as a rejection probability 
P r . Another aspect we monitor is the frequency of generating configurations with nearly 
singular matrices Mj(s, L t a t ). For this we introduce a small positive parameter S and reject 
configurations with 



|detMj(s,L t oi t )| < 5 



N 



i=l,...,N 



(74) 



By taking the limit 5 — > + we can determine if poorly-conditioned matrices make any 
detectable contribution to the final result. The rejection rate due to these singular matrices 
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is included in the total rejection probability P T , but is also recorded separately as a singular 
matrix probability P s . 

We use double precision arithmetic in all calculations. Numerical stabilization methods 
based on Gram-Schmidt orthogonalization have been developed for determinantal quantum 
Monte Carlo simulations 55, 0, [57)]. These methods have proved successful in reducing 
round-off error for singular matrix determinants and should also be effective here. However 
in our ground state lattice simulations we find that round-off error is only one of several 
problems associated with singular matrix configurations. Another problem that arises is 
quasi-non-ergodic behavior in the Monte Carlo updates. For the case of hybrid Monte 
Carlo importance sampling this problem can be visualized in terms of the landscape of the 
function Hj(s,p) defined in Eq. (j5T?|) . Each singular matrix configuration s produces a sharp 
peak in Hj(s,p). Since hybrid Monte Carlo trajectories follow the contour lines of Hj(s,p), 
these trajectories can get trapped on orbits around singular matrix configurations. Since 
the matrix inverse of Mj diverges near these singular configurations, this behavior is also 
accompanied by large fluctuations in the difermion spatial correlation function 

(a\ (n)a{(n)a T (0H(0)}. (75) 



This problem was observed in lattice simulations in [42|]. For these reasons we take a 
cautious approach to singular matrix configurations until the problem is better understood. 
We use the singular matrix probability P s as a diagnostic tool to measure the severity of the 
problem. 

We use the unpolarized ten-particle system to test the computational performance of 
the four auxiliary-field methods. We compute £5,5(^0%) for spatial length L = 5 and 
temporal lengths L t = 24, 48, 72. The simulations are performed using 8 processors starting 
with different random number seeds producing independent configurations. As before the 
individual processor results are averaged, and the standard deviation is used to estimate 
the error of the average. For j = 1,2,4 we generate hybrid Monte Carlo trajectories with 
N s te P = 10 and e step = 0.1. For the j = 3 Metropolis update we flip the sign of s for 0.15% 
of lattice sites selected randomly on each lattice sweep. For the singular matrix condition 
in Eq. flHD, we use 5 = 5 x 10~ 7 . 

Results for ^ 5t5 (L t a t ) for L t = 24 are shown in Fig. [31 For j = 1,2,4 the horizontal 
axis is the number of completed hybrid Monte Carlo trajectories per processor. For j = 3 
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FIG. 3: (Color online) Monte Carlo results for £5 s(Ltat) for L = 5 and Lt = 24. P r is the rejection 
probability, and P s is the singular matrix probability. 



the horizontal axis is the number of attempted Metropolis updates per lattice site. On an 
Intel Xeon processor the CPU run time for 1000 hybrid Monte Carlo trajectories is roughly 
the same as 150,000 lattice sweeps with Metropolis updating on 0.15% of the sites. This 
corresponds with 225 updates per site. 

All four calculations are in agreement with an average value £ 5j5 (24a t ) m 0.39. In all 
cases the estimated errors are relatively small. In order of increasing error, the error bars 
for j = 4 is smallest, then j = 2, j = 1, and j = 3. This is the same ordering we get by 
sorting rejection probability P r from lowest to highest. The singular matrix probability P s 
is zero in all cases. 

Results for L t = 48 are shown in Fig. HI This time all four calculations are in agreement 
with £5,5(480^) ~ 0.29. In all cases the estimated errors are still rather small, but upon 
closer inspection there are some early signs of trouble for j = 3 and j = 2. For j '• = 3 
the rejection probability is quite large. This suggests difficulties in sampling the space 
of discrete auxiliary-field configurations. For j = 2 the singular matrix probability P s is 
1.9%. This is approaching the level where the contribution due to singular matrices may 
be detectable above the background of stochastic error. 

Results for L t = 72 are shown in Fig. [5j The simulation for j = 2 has failed due to the 
high singular matrix probability. The calculations for j = 1,3,4 are consistent with a value 
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FIG. 4: (Color online) Monte Carlo results for £5,5(^0^) for L = 5 and L t = 48. 



0.4 
0.3 - 
0.2 [1 
0.1 




7=1 
n — 1 — 1 — r 



J I I I L 



2 4 6 8 10 

1000 HMC traj. 
P r = 13.9% 
P s = 0.2% 

(a) 



0.4 



0.3 -x-.-... 



0.2 \* 
0.1 




7 = 2 
n — 1 — 1 — r 



1 I i"i : 1 



2 4 6 8 10 

1000 HMC traj. 
P r = 19.9% 
P s = 13.0% 
(b) 



0.4 
0.3 



0.1 




7 = 3 
n — 1 — 1 — r 



0.2 H* '"^%!HS*^- 



J I I I L 



2 4 6 8 10 



0.3 - 




2 4 6 8 10 



225 Met. updates/site 1000 HMC traj. 



P r = 37.6% 
P s = 0.0% 



P r = 5.1% 
P s = 0.0% 



(c) 



(d) 



FIG. 5: (Color online) Monte Carlo results for ^^{L t at) for L = 5 and L t = 72. 



of £5,5(72*^) pa 0.25. However the rejection probability for j = 3 is very high. This could 
explain the gradual downward and then upward drift in the j = 3 data, a sign that the 
Monte Carlo simulation has not fully equilibrated due to a long autocorrelation time. The 
same is true to a lesser extent for j = 1. For j = 1 we also see that P s is small but nonzero. 
This suggests that problems with singular matrices may appear for somewhat larger L and 
L t . 
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TABLE III: Lattice dimensions L 3 x Lt used in calculations for £5,5 (-^t Q t)- 



T 3 


4 3 


c:3 

(J 


p3 
L) 


7 3 


8 3 




16 


24 


36 


48 


70 


u 


20 


30 


42 


56 


80 




48 


72 


102 


144 


190 



TABLE IV: Lattice dimensions L 3 x Lt used in calculations for ^Tj(Ltat). 



L 3 


4 3 


5 3 


6 3 


7 3 


8 3 




12 


18 


30 


40 


50 


it 


16 


24 


36 


48 


80 




44 


72 


96 


128 


170 



From this analysis we rate the performance for j '• = 4 superior to the other three methods. 
This is not entirely unexpected. The use of a bounded auxiliary field should reduce the 
likelihood of exceptional configurations producing singular matrices. The use of a continuous 
auxiliary field makes it possible to use hybrid Monte Carlo which is better at reducing 
rejection probability than local Metropolis updates. For fixed trajectory length iV stC p£: stC p, 
the rejection probability for hybrid Monte Carlo scales quadratically with step size, e^ tep . 
This is in contrast with the Metropolis algorithm, where the rejection probability scales 
linearly with the fraction of lattice sites updated in each sweep. This contributes to a much 
slower performance of the Metropolis algorithm at large L and L t . 

VI. MAIN RESULTS 

We use the bounded continuous auxiliary-field formulation j = 4 for the lattice results 
presented in this section. We consider the unpolarized ten-particle and the fourteen-particle 
systems. For the calculation of £ 5j5 (L 4 a 4 ) we use the lattice dimensions L 3 x L t shown in 
Table ED For £ 7J (L t a t ) we use the lattice dimensions L 3 x L t in Table IIVI We use 
trajectory parameters N step = 10, e step = 0.1, and singular matrix parameter 5 = 5 x 10 -7 . 
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The rejection probability P r reaches a maximum of 12% and the singular matrix probability 
P s reaches a maximum of 3% for the largest lattice volume, L 3 = 8 3 and L t > 170. For 
most of the lattice simulations the values for P T and P s are much smaller. The simulations 
are run with a minimum of 16 processors each running a minimum of 10, 000 hybrid Monte 
Carlo trajectories. 



From Eq. ( 1531) we expect 

UnH) ~ Un + ( E ^- E ° N A e-K-*V> (76) 

k^o\ C N,N\ \ ^*,N J 

at large t. To extract £n,n we perform a least squares fit of £jv,jv(£) to the functional form, 

UN(t) = Un + b NtN e~ s ^ Ept . (77) 

For asymptotically large E F t, we can identify 5n,nE f with the energy separation between 
the ground state and the first excited state, possibly with degenerate partners. For large 
N the lowest excitation is expected to be a two phonon state with zero total momentum. 
For large N the excitation energy for this state is small compared with Ep and therefore 
Sn,n « 1- 

Our observable £j\r,iv(X) is proportional to the expectation value of the energy. When the 
energy difference N — E^ N is very small compared to N the contribution proportional 
to e - ^^'^ - ^'^) 4 is difficult to resolve against the background of stochastic noise. We note 
the factor 

N,N ^N^N /-oN 

E~° ( } 

111 N,N 

multiplying e - ^^^ - ^,^)* j n Eq. fl76|) . If the objective were to measure N — E^ N 
accurately for very low excitations, then it would be more effective to compute the Euclidean 
time projection of some other observable such as the difermion spatial correlation function 

(a|(n)a{(n)a T (0) ai (0)). (79) 

This technique was used to identify the lowest energy excitations for several unpolarized 
lattice systems at unitarity j^ . 

In the analysis here we focus only on measuring the ground state energy accurately and 
ignore the numerically small contributions hidden in the asymptotic tail of ^jvat(^)- We 
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TABLE V: Results for the three-parameter fit for £5,5(2) using 65,5, £5,5, £5,5. 



L 


U D,0 


u 0.0 




y 2 /d f 


4 


0.42(12) 


0.47(13) 


0.22(2) 


1.2 


5 


0.40(7) 


0.44(7) 


0.23(1) 


0.4 


6 


0.41(8) 


0.42(8) 


0.24(1) 


1.6 


7 


0.81(18) 


0.69(10) 


0.26(1) 


0.6 


8 


0.36(12) 


0.32(15) 


0.24(4) 


0.6 



TABLE VI: Results for the three-parameter fit for £7,7(2) using 67,7, £7,7, £7,7. 



L 


£>7,7 


£7,7 


£7,7 


X 2 /d± 


4 


0.24(9) 


0.39(13) 


0.26(1) 


1.3 


5 


0.28(4) 


0.35(6) 


0.27(1) 


1.0 


6 


0.39(7) 


0.45(6) 


0.29(1) 


1.3 


7 


0.29(8) 


0.27(15) 


0.26(4) 


1.7 


8 


0.36(7) 


0.41(10) 


0.30(2) 


0.8 



determine &/v,iv> $n,n, £jv,jv from least squares fitting over the range Ept = 2 to Ept = 9. 
The values 6jv,at and <5jv,v we determine from least squares fitting should be interpreted as 
a spectral average, 

£ j^rf ( <,iV -n,AA g-^-^), w 6 ^ e -«zr^. (80) 

k^0 \c° N ,n\ \ E N,N J 

Since £,N,N{t) has a well-defined continuum limit for fixed Ept, each of the dimensionless 
parameters 6jv,at, <^v,at, £a?,v also has a well-defined continuum limit. Up to uncertainties the 
size of least squares fitting errors, £n,n is independent of the initial state overlap amplitudes 
c k NN . However b N N and both depend on c k N ^ N . Table IVl shows the three-parameter 
fit results for N = 5, and Table [VTl shows the three-parameter fit results for N = 7. The 
average chi-square per degree of freedom for the fits is about 1. The error estimates for 
the fit parameters are calculated by explicit simulation. We introduce Gaussian-random 
noise scaled by the error bars of each data point for £at,7v(2). The fit is repeated many 
times with the random noise included to estimate the one standard-deviation spread in the 
fit parameters. 
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TABLE VII: Results for the two-parameter fit for £5,5^) using 65,5, £5,5 with £5,5 = 0.47. 



L 




r* EC c; 
SO,0 


Y 2 /d f 


4 


0.42(4) 


0.223(5) 


1.0 


5 


0.43(4) 


0.236(2) 


0.3 


6 


0.47(3) 


0.247(2) 


1.5 


7 


0.51(4) 


0.242(5) 


0.9 


8 


0.47(4) 


0.264(6) 


0.6 



TABLE VIII: Results for the two-parameter fit for £7,7(0 using 67,7, £7,7 with £7^ = 0.37. 



L 


^7,7 


£7,7 


rVd.f. 


4 


0.23(2) 


0.261(4) 


1.1 


5 


0.29(1) 


0.276(3) 


0.9 


6 


0.31(2) 


0.285(2) 


1.4 


7 


0.34(4) 


0.282(9) 


1.6 


8 


0.33(3) 


0.294(6) 


0.8 



The error in £n,n would be considerably smaller if the fit needed only two parameters 
rather than three parameters. This can be arranged if we neglect the L-dependence of 
5n,n and fix 5n,n according to the average values for L — 4,5, 6, 7, 8 as quoted in Tables IVl 
and [VI] Since 5n,n has a well-defined continuum limit, the L-dependence of Sn,n should 
in fact be small. For N = 5 the average value is $5,5 = 0.47, and for N = 7 the average 
value is 5j <7 = 0.37. Using these values we refit £jv,jv (t) with the two parameters 6jv,AT) 6v,at- 
Table IVIII shows the two-parameter fit results for N = 5, and Table IVIIII shows the two- 
parameter fit results for N = 7. 

The average chi-square per degree of freedom is again about 1. The lattice data for 
£5,5^) together with the two-parameter fit functions are shown Fig. [61 Fig. [7] shows the 
lattice data for ^7,7 (t) with two-parameter fit functions. The lattice results show that £5,5(0 
and £7,7(0 are both approximately universal functions of Ept independent of L. This was 
expected from the scale invariance of the unitarity limit. 

Using the results for £ 5i5 and £7,7 for L = 4,5, 6, 7, 8 we can extrapolate to the continuum 
limit L — > 00. We expect some residual dependence on L proportional to L~ x , arising from 
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FIG. 6: (Color online) The lattice data for £5,5(2) versus for L = 4,5,6,7,8. Also shown are 
the results of the two-parameter fits. 



effects such as the effective range correction, broken Galilean invariance, and possibly other 
lattice cutoff effects. From the three-parameter fit results in Tables |V] and I VI I the linear 
extrapolation in L _1 gives the continuum limit values 

£5,5 = 0.308(25), (81) 

£ 7J = 0.337(30). (82) 
If we extrapolate the two-parameter fit results in Tables IVHI and IVIIII we get 

£5,5 = 0.292(12), (83) 

£ 7>7 = 0.329(5). (84) 

Results from the two-parameter fits for £ 5 5 and £7,7 at finite L and the corresponding contin- 
uum limit extrapolations are shown in Fig. [HI We note that the continuum limit fit for ^n,n 
could also be performed with E ^^ defined in the continuum limit for the same cubic box 
size. This procedure is not recommended since it introduces a larger L~ 2 dependence and 
degrades the quality of the linear L" 1 fit. However the fit can be done and the extrapolated 
values for £ 5i5 and £7,7 are each about 0.015 higher than the values reported in Eq. (183!) and 
fl84"j) with somewhat larger error bars. 
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FIG. 7: (Color online) The lattice data for £7,7^) versus Ept for L = 4,5,6,7,8. Also shown are 
the results of the two-parameter fits. 
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FIG. 8: (Color online) Results from the two-parameter fits for £5^ and £7^ at finite L and the 
corresponding continuum limit extrapolations. 



VII. DISCUSSION 



In 40| the values 



£5,5 = 0.24(2), 



(85) 
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6,7 = 0.28(2), (86) 

were found based on an average of lattice results for L = 4,5, 6. The results here for £5^ and 
£ 7 7 with L — 4, 5, 6 are consistent with these values. The continuum limit extrapolation was 
not possible in [40j due to problems with increasing singular matrix probability P s . That 
calculation used the auxiliary field formulation j = 2, which had the largest P s of the four 
methods considered here. The bounded continuous auxiliary-field method appears to solve 
this problem for the lattice systems considered here. The values for £ 5 5 and £7,7 from small 
lattices L = 4, 5, 6 are each shifted upwards by 0.05 when extrapolated to the continuum 
limit. 

In the calculations presented here we have only considered systems with 10 and 14 par- 
ticles and have not attempted to determine the thermodynamical limit iV — > 00. It is 
therefore interesting to compare with results obtained using other methods that have com- 
puted ground state energies for both small and large values of N. Each of these other 
methods contain some unknown systematic errors, and so a benchmark comparison with 
continuum extrapolated Monte Carlo lattice results for N = 5 and N = 7 provides an 
estimate of the systematic error. 

There is a discrepancy of about 0.13 in the reported values for £ 5 5 and £7,7 between 
continuum extrapolated lattice results and fixed-node Green's function Monte Carlo results 
on a periodic cube. The fixed-node Green's function Monte Carlo simulations find £n,n — 



0.44(1) for 5 < iV < 21 [27|] and 0.42(1) for larger N |28|, |23|. This discrepancy suggests 
that the upper bound on the ground state energy using fixed-node Green's function Monte 
Carlo might be lowered further by a more optimal fermionic nodal surface. 

A recent Hamilton.au lattieeatndy eompnted ^ for the range TV — 2 to TV — 32 on 
cubic lattices up to L 3 = 16 3 [4l|]. This calculation used a method called the symmetric 
heavy-light ansatz in the lowest filling approximation. In Fig. [9] we compare Monte Carlo 
lattice results presented here and symmetric heavy-light ansatz results for N = 5 and N = 7. 
The continuum limit extrapolations of the symmetric heavy-light results give 

£ 5 , 5 = 0.301(1), (87) 

£ 7 , 7 = 0.339(1). (88) 
These results are within 0.01 of the values found in Eq. ( 1831) and ( 1841) . This level of accuracy 
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FIG. 9: (Color online) Comparison of results from Fig. [8] with Hamiltonian lattice results using 
the symmetric heavy-light ansatz in the lowest filling approximation. 



is consistent with the size of errors found in 



411 ] for four-body and six-body systems of the 



ID, 2D, 3D attractive Hubbard models at arbitrary coupling. 

The different L~ l slopes for the Hamiltonian lattice and Euclidean lattice extrapolations 
in Fig. [9] are consistent with the fact that the effective range for the Hamiltonian lattice 
interaction is a smaller negative fraction of the lattice spacing than for the Euclidean lattice 
transfer matrix. However there are other effects such as broken Galilean invariance on the 
lattice which produce a similar L _1 dependence 



. An accurate calculation of the effective 



range correction with controlled systematic errors requires either an effective range larger 
than the lattice spacing or an analysis of the effect of changing the effective range parameter 
relative to the lattice spacing. The effective range correction has recently been computed 
on the lattice with realistic dilute neutron matter at next-to-leading order in chiral effective 
field theory |5a,lfiQ|. 

Results from the symmetric heavy-light ansatz for general N are shown in Fig. [ID] j^l]. 
We note that the value for £n,n reaches a maximum at iV = 7. This can be explained by 
the closed shell at iV = 7 in the free fermion ground state and the absence of shell effects 
in the interacting system. The numerical agreement for the benchmark comparisons at 
N = 5 and N = 7 provides some confidence in the symmetric heavy-light ansatz value of 
£ = 0.31(1) for the unitarity limit in the continuum and thermodynamic limits [41]. 
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The bounded continuous auxiliary-field method succeeds in reducing the singular ma- 
trix probability P s by minimizing large fluctuations in the auxiliary-field transfer matrix 
elements. Very roughly this corresponds with reducing the size of 



max \Aj [s (n,n t )} \ , (89) 

s{ft,nt) 



given the constraints 

Jd jS {n,n t ) 1 = 1, (90) 
djS (ft, n t ) Aj [s (n, n t )} = 0, (91) 

djS (ft, n t ) A 2 j [s (n, n t )\ = -Ca t . (92) 

For lattice systems larger than the ones considered here the problem with singular matrices 
will reappear. For very large systems there may be no choice but to use the stabilization 



methods developed in |55|, |56|, |57| to reduce round-off error and confront the problems of 
large fluctuations and quasi-non-ergodic behavior using brute-force large-scale simulations. 
However for moderately larger systems it may be possible to gain further advantage using 
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a function Aj [s(n, n t )] with a steep slope at s(n,n t ) = and relatively flat away from 
zero. The most extreme case would be to use an odd periodic step function. But this 
is equivalent to the discrete auxiliary-field formulation j — 3 with the problems discussed 
concerning large rejection probabilities. However with a smooth approximation to an odd 
periodic step function, it may be possible to reduce the singular matrix probability while 
at the same time compensate for the increase in rejection probability with a smaller hybrid 
Monte Carlo step size, e stC p- 

VIII. SUMMARY 

We have presented new Euclidean lattice methods which remove some computational bar- 
riers encountered in previous lattice calculations of the ground state energy in the unitarity 
limit. We compared the performance of four different auxiliary-field methods that produce 
exactly the same lattice transfer matrix. By far the best performance was obtained using a 
bounded continuous auxiliary field with hybrid Monte Carlo updating. With this method 
we calculated results for 10 and 14 fermions at lattice volumes 4 3 , 5 3 , 6 3 , 7 3 , 8 3 and extrapo- 
lated to the continuum limit. For 10 fermions in a periodic cube, we found the ground state 
energy to be 0.292(12) times the ground state energy for non-interacting fermions. For 14 
fermions the ratio is 0.329(5). These values may be useful as benchmarks for calculations 
of the unitarity limit ground state using other methods. 
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